function rho = newton(nd,rho,omsc,gamsc,z,cset,coef) 
 % L rho = raux1

      raux1 =  liouv(nd,rho,omsc,gamsc,cset,coef);

% raux1 = L rho - z 1op rho.

            raux1=raux1-z*rho;
            rho=raux1;

      end
